


%% 信道系数的优化
function SignalPower_xbyb = NearFieldFDA2D(x_h, y_h, D, lambda)
% 定义参数和变量
N_d = floor(2 * D / lambda); % 垂直方向的阵元数  600
N_e = floor(2 * D / lambda); % 水平方向的阵元数  600
spacing = D / N_d;
M = 2; % 用户数
b = 2; % Boresight gain参数
% x_h = 5; % 用户 Bob 的 x 坐标
% y_h = 5; % 用户 Bob 的 y 坐标
z_h = 0; % 用户 Bob 的 z 坐标

k = 2 * pi / lambda; % 波数
alpha_c = 1.2;
beta_c = 827.67;



% 生成信道系数
% rho = rand(N_d, N_e); % 元件位置
A_h = zeros(N_d, N_e);  % 这个矩阵没有意义，只是暂存数据,暂存各个天线算出来的数据
A_b = zeros(N_d, N_e);

s = randn(N_d, N_e); % 模拟发送信号

SignalPower_xbyb = 0;


% p_i_l = [(i - 1) * spacing - D / 2, 0, (l - 1) * spacing - D / 2];
% p_b = [x_b, y_b, z_b];
% Theta_i_l_b = calculateTheta(p_i_l, p_b);
% distance_b_il = norm(p_b - p_i_l); % 计算距离 |p_m - p_i_l|
% phase_b_il = exp(-1i * k * distance_b_il); % 计算相位 e^{-j*k*|p_m - p_i_l|}
% A_b(i, l) = sqrt(F(Theta_i_l_b)) * (lambda / (4 * pi * norm(p_b - p_i_l))) * phase_b_il;


for i = 1:N_d
    for l = 1:N_e
        % rho(i, l) = 
        p_i_l = [(i - 1) * spacing - D / 2, 0, (l - 1) * spacing - D / 2];
        % p_i_l = [0, (i - 1) * spacing - D / 2, (l - 1) * spacing - D / 2];
        % 计算用户 Bob 的信道系数
        p_h = [x_h, y_h, z_h];
        Theta_i_l_h = calculateTheta(p_i_l, p_h);
        distance_h_il = norm(p_h - p_i_l); % 计算距离 |p_m - p_i_l|
        phase_h_il = exp(-1i * k * distance_h_il); % 计算相位 e^{-j*k*|p_m - p_i_l|}
        A_h(i, l) = sqrt(F(Theta_i_l_h)) * (lambda / (4 * pi * norm(p_h - p_i_l))) * phase_h_il;
        % SignalPower_xbyb = SignalPower_xbyb + A_h(i, l);
    end
end

A_vec = A_h(:);  % 列向量，Nt * Nr

[U, S, V] = svd(A_vec);
% 选择前 N_s 个奇异值，这里假设 N_s = 1
N_s = 1;
U_s = U(:, 1:N_s);
V_s = V(:, 1:N_s);

a_beam = A_vec' * V_s;   % Nr * Nt * Nt * d;  % Nr * d;
% 应该要求的H为 Nr * Nt
% SignalPower_xbyb = norm(SignalPower_xbyb)^2;
SignalPower_xbyb = norm(a_beam)^2;

end
% 

% 
% %% DMA的发送信号
% N = N_d * N_e; % 总元件数量
% sigma_b2 = 1; % 噪声方差
% 
% % 生成随机的数字前向编码向量
% w_b = (randn(N_d, 1) + 1i * randn(N_d, 1));
% w_w = (randn(N_d, 1) + 1i * randn(N_d, 1));
% 
% 
% % 初始化 H 矩阵和 Q 矩阵
% H = zeros(N, N);
% Q = zeros(N, N_d);
% 
% % 初始化示例的元件相位和位置
% 
% 
% phi = linspace(0, 2 * pi, N_d); % 元件响应角   待优化的变量
% 
% % 填充 H 矩阵和 Q 矩阵
% for i = 1:N_d
%     for l = 1:N_e
%         % 计算元件相位
%         
%         h_i_l = exp(  ((l-1)*spacing-D/2)  *  (alpha_c + 1i*beta_c)  );
%         
%         % 计算可调响应 q_i_l
%         q_i_l = (j + exp(1i * phi(i))) / 2;
%         
%         % 将 h_i_l 填充到 H 矩阵对角线
%         H((i - 1) * N_e + l, (i - 1) * N_e + l) = h_i_l;
%         
%         % 将 q_i_l 填充到 Q 矩阵
%         Q((i - 1) * N_e + l, i) = q_i_l;
%     end
% end
% 
% % 生成示例输入信号 x_h
% x_h = ( randn(1, 1) + 1i * randn(1, 1) ) / sqrt(2);
% x_w = ( randn(1, 1) + 1i * randn(1, 1) ) / sqrt(2);
% 
% s = H * Q * (w_b * x_h + w_w * x_w);
% 
% % % 计算可达速率
% % a_m = ones(1, N); % 示例 a_m 向量
% % w_j = ones(N, 1); % 示例 w_j 向量
% % R_m = log2(1 + abs(a_m * (H * Q * w_j))^2 / (sum(abs(a_m * (H * Q * w_j))^2) + sigma2));
% % 
% % disp(['可达速率 R_m = ' num2str(R_m)]);
% 
% s = s(:); % 将发送信号转化为列向量
% n_b = sqrt(0.5) * (randn(1, 1) + 1i * randn(1, 1)); % 用户 Bob 的噪声
% n_w = sqrt(0.5) * (randn(1, 1) + 1i * randn(1, 1)); % 窃听者 Willie 的噪声
% r_b = a_b' * s + n_b; % 用户 Bob 的接收信号  
% r_w = a_w' * s + n_w; % 窃听者 Willie 的接收信号
% 
% disp('用户 Bob 的接收信号:');
% disp(r_b);
% 
% disp('窃听者 Willie 的接收信号:');
% disp(r_w);
% 
% 
% %% 待优化变量： Pb Q w_b w_w
% Pbmax = 10;
% Pb = 1; % 只剩Pb没初始化了， variable Bob的发射功率待优化
% epsilon = 0.01;
% 
% % 生成一个在[0, Pwmax]范围内的均匀分布的Pw
% Pwmax = 2;
% Pw = Pwmax * rand;
% 
% 
% % Bob处接收到的信噪比
% gamma_b = Pb * norm(a_b' * H * Q * w_b)^2 / ( Pw * norm(a_b' * H * Q * w_w)^2 +  sigma_b2) ;  % aim to max this
% disp('Bob 的信噪比:');
% disp(gamma_b);
% 
% gamma_w = Pwmax/2 * norm(a_w' * H * Q * w_w)^2 / ( Pw * norm(a_w' * H * Q * w_w)^2 +  sigma_b2) ;  % aim to max this
% disp('Bob 的信噪比:');
% disp(gamma_b);
% 
% % 
% % 
% % % constrains
% % % Decision Error Probility of Willie
% % Pr_b = 1 - Pb * norm(a_w' * H * Q * w_b)^2 / ( Pwmax * norm(a_w' * H * Q * w_w)^2 );
% % 
% % Pr_b >=  epsilon;
% % 
% % Pb <= Pbmax;
% % 
% % norm(w_b) == 1;
% % 
% % norm(w_w) == 1;
% % 
% % % Q的约束
% % 
% 
% 
% 
% 
% 
% 



